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A  FORTRAN  SUBROUTINE  TO  EVALUATE 
SPHERICAL  BESSEL  FUNCTIONS  OF  THE  FIRST,  SECOND, 
AND  THIRD  KINDS  FOR  COMPLEX  ARGUMENTS 


1.  INTRODUCTION 

^  To  determine  analytically  the  acoustic  scattering  from  absorbing  spheres,  it  is  necessary  to  evalu¬ 
ate  spherical  Bessel  functions  of  the  first,  second,  and  third  kinds  for  complex  arguments,  covering  a 
range  of  integer  orders  from  zero  to  a  value  approximately  equal  to  the  absolute  value  of  the  argument. 
A  double-precision  (G-format)  FORTRAN  subroutine,  CSPJYD,  has  been  written  for  the  VAX-750 
that  will  generate  tables^of  jn(z)  and  y„(z)  and/or  /j„u'(z)  for  a  given  value  of  z  and  a  range  of  n  from 
zero  to  a  given  maximum  it."  . 

2.  METHODS  OF  COMPUTATION 

For  an  absolute  value  of  the  argument  less  than  or  equal  to  0.5,  j„(z)  and  y„(z)  are  each  gen¬ 
erated  by  an  alternating  series;  they  are  then  used  to  generate  the  //BU)(z)  values.  For  arguments  of 
absolute  value  greater  than  0.5,  when  the  absolute  value  of  the  coefficient  of  the  imaginary  part  of  the 
argument  is  less  than  5.0,  j„(z)  and  y„(z)  are  always  calculated,  and  conditional  on  the  subroutine  call, 
only  then  is  /;n(*'(z)  derived.  However,  when  the  absolute  value  of  the  coefficient  of  the  imaginary  part 
of  the  argument  is  greater  than  or  equal  to  5.0,  j„(z)  and  /t„u)(z)  are  always  calculated,  and  yB(z)  only 
if  requested  in  the  subroutine  call.  The  values  for  y„(z ),  or  /j„u)(z),  will  not  be  returned  to  the  calling 
program  unless  the  user  so  requests,  even  though  those  functions  were  calculated.  See  Sec.  4,  Instruc¬ 
tions  for  Use,  for  further  clarification. 


In  the  following  equations,  z  -  u  +  /v  was  used  in  place  of  the  more  common  notation 
z  -  v  +  />,  to  prevent  confusing  the  coefficient  of  the  imaginary  part  of  the  argument  with  the  spheri¬ 
cal  Bessel  function  of  the  second  kind,  >„(z). 


a. 


if  I z t  0.5, 


13-5 
—  13  5 


(2 n  +  1) 

•  •  (2  n-  1) 


(z2/  2) 


1  - 


l!<2  n  +  3) 
(z2/  2) 


y"(:)  **'  '*  I !  ( l  —  2n) 

/t,'n(z)  —  j„(z)  +  (y„(z)j  /»„M)(z)  is  calculated  if  v 
/t„,2l(z)  —  j„(z )  -  t>„(z)|  i/„<2)(z)  is  calculated  if  v 


(z2/  2)2 

2!(2n  +  3)  (2 n  +  5) 

,  (z72)2 

2KI  -  2 n)  (3  -  2 n) 

£  0 
<  0 


b. 


if  I z I  >  0.5  and  |v|  <  5.0, 
sin  z 


jo  <z) 
7i(z)  ■ 
To<z)  - 
y,(z)  - 


sin  z  _  cos  z 
Z2  z 

-  cos  z 


-  cos  z 


sin  z 


where: 

sin  z  -  sin  w(ev  +  e~')/2  +  i  (cos  u  (ev  -  e~v)/2]  and 
cos  z  -  cos  «(ev  +  e~')/2  -  /Isin  u(e'  -  e~')/2] 
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yB(z)  “  (2/i  +  3 )z~l  yB+1(z)  -  jH+2(z)  (backward  recursion) 
y„(z)  —  (2n  -  l)z-1  y„-i(z)  —  y„~2 (z)  (forward  recursion) 
/i„<n (2)  -  j„(z)  +  /yB(z)l  hj,}(z)  is  calculated  if  v  ^  0 

A„<2,(z)  -  j„(z)  -  iy„(z) J  h „(2)(z)  is  calculated  if  v  <  0 


c.  if  \z |  >  0.5  and  } v)  >  5.0, 
sin  z 


z 

sin  z 


cos  z 


where  sin  z  and  cos  z  are 
defined  as  in  b. 


jo(z) 

'|(2)  z' 

/»o(1)  (z)  ”  e_v(sin  w  -  1  cos  «)/z 

/ti(1)  (z)  ■■  e~v(sin  u  —  /  cos  u)/z2—  e~y  (cos  u  +  /  sin  u)/z| 
/r0(2)  (z)  -  ev(sin  «  +  i  cos  w)/z 
/j/2*  (z)  —  ev(sin  u  +  /  cos  u)/z2—  ev(cos  u  -  i  sin  u)/z 


hk"  (z)  and  A(<u  (z)  are 
calculated  if  v  >  0 


/»o<2)  (z)  and  /i1(2)  (z)  are 
calculated  if  v  <  0 
yB(z)  -  (2/i  +  3)z-1  yB+1(z)  -  j„+2(z)  (backward  recursion) 

ABa>(z)  -  (2/i  -  l)z-1  /»„**)  (z)  —  h^\  (z)  (forward  recursion) 

.VB(z)  -  (/B(z)  -  //tB0)(z)  or 
y,(z)  -  -  //„(z)  +  ih}2Hz) 


3.  VERIFICATION 


Because  published  tables  often  are  limited  in  range  or  lack  precision,  the  accuracy  of  the  derived 
function  values  was  checked  with  at  least  one  of  the  following  Wronskian  relationships: 

a-  A+t(z)  >B(z)  -  7„(z)  JV*. i(z)  -  1/z2 

b.  U(z)  //ji’i  (z)  -  7b+i(z)  /tB(,)(z)]  /  -  1/z2 

c.  0B (z)  b„(}\  (z)  -  yB+1(z)  /tB<2,(z)J/  -  1/z2 


The  table  below  lists  some  of  the  arguments  and  orders  for  which  runs  were  made,  the  degree  of 
accuracy  in  the  resultant  Bessel  functions,*  and  their  respective  Wronskians.  Note  that  these  Wronski- 
ans  were  calculated  as  functions  of  j„-\(z),  jn(z),  h^t\  (z),  and  hjk)(z' ,  in  every  case. 


Argument 

max  n 

Accuracy  of  Results 
at  max  n  (places) 

Wronskian 
Agreement  at  max  n 
(figures;  places) 

0.01 

— 

O.OOli 

20 

10 

14;  10 

1.0 

- 

lOO.Oi 

100 

10 

1 4;  1 9 

100.0 

± 

0.5i 

100 

10 

13;  1 7 

100.0 

- 

lOO.Oi 

100 

10 

16;20 

1000.0 

- 

lO.Oi 

100 

10 

11,18 

1000.0 

— 

lOO.Oi 

100 

10 

15;21 

0.0 

± 

0.4i 

10 

(a) 

16;  1 5 

0.0 

± 

0.6i 

10 

(a) 

1 5;  1 4 

0.0 

± 

5.1i 

10 

(a) 

15;16 

’’’No  comparison  was  made. 


•The  accuracy  was  determined  by  comparing  these  results  with  a  10-place  table  of  spherical  Bessel  functions  of  the  first  and 
second  kinds. 
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4.  INSTRUCTIONS  FOR  USE 

The  subroutine  call  is: 

CALL  CSPJYD  (AR,  AI,  N,  SJR,  SJI,  SYR,  SYI,  SHR,  SHI,  NAK) 

input:  AR  —  the  real  part  of  the  argument  z 

AI  —  the  imaginary  part  of  the  argument  z 
N  “  the  maximum  order. 

NAK  ~  2,  if  j'(z)  and  y„(z )  are  to  be  calculated. 

NAK  “  3,  if  j„(z)  and  h}k)(z)  are  to  be  calculated. 

NAK  —  S,  if  all  three  are  to  be  calculated. 

output:  SJR  —  a  table,  from  n  -  0  to  n  -»  N,  of  the  real  part  of  J„(z)- 

SJI  —  a  table,  from  n  —  0  to  n  *  N,  of  the  imaginary  part  of  jH(z). 

SYR  —  a  table,  from  —  0  to  n  —  V,  of  the  real  part  of  y„  (z)  if  NAK  ™  2  or  5. 

SYR  —  a  table,  from  n  —  0  to  n  »  N,  of  the  real  part  of  h%k)(z)  if  NAK  ■»  3. 

SYI  —  a  table,  from  n  0  to  /r  »  N,  of  the  imaginary  part  of  y„(z)  if  NAK  —  2  or  5. 

SYI  -*  a  table,  from  n  -  0  to  n  -  N,  of  the  imaginary  part  of  h}kHz)  if  NAK  “  3. 

SHR  —  a  table,  from  n  -  0  to  n  «■  N,  of  the  real  part  of  h„(k>(z)  if  NAK  *■  5. 

SHI  <-  a  table,  from  n  -  0  to  n  »  N,  of  the  imaginary  part  of  h„(k>(z)  if  NAK  —  5. 

If  NAK  —  2  or  3,  for  SHR  and  SHI  use  dummy  parameters,  which  do  not  have  to  be  dimen¬ 
sioned.  All  parameters  except  N  and  NAK  must  be  REAL*8.  The  routine  calculates  (u  +  »'v)  for 
positive  v  and  / ij2)(u  +  (v)  for  negative  v.  SJR  and  SJI  should  be  dimensioned  by  at  least 
( u 2  +  v2)1/J  +  30  or  N  +  30,  whichever  is  the  larger.  SYR  and  SYI  should  be  dimensioned  by  at  least 
N  +  1.  If  NAK  —  2  or  3,  SHR  and  SHI  require  no  dimensioning,  but  if  NAK  —  S,  they  also  should 
be  dimensioned  by  at  least  N  +  1. 

Subroutine  CSPJYD  calls  two  other  subroutines,  DVDD  and  MLTD;  they  perform  complex  divi¬ 
sion  and  multiplication,  respectively. 

5.  PORTABILITY 

Generally  speaking,  CSPJYD  can  be  adapted  for  use  on  many  other  computers.  It  works  with 
integers  and  double-precision  real  values  only,  so  there  is  no  problem  in  using  it  on  a  computer,  such 
as  the  PDP-11,  which  permits  no  higher  precision  for  its  complex  numbers  than  COMPLEX*8.  The 
checks  for  •  srflow,  divide  by  zero,  and  underflow,  which  are  located  in  subroutines  CSPJYD  and 
DVDD,  would  have  to  be  changed  for  non-DEC  machines,  but  this  should  not  prove  difficult. 

6.  LISINGS  AND  OUTPUT1 

Following  are  source  listings  of  subroutines  CSPJYD,  DVDD,  and  MLTD;  a  listing  of  the  test 
program  TSPHBF;  and  output  from  two  sample  runs  of  program  TSPHBF.  One  note  here:  if  both 
y,(z)  and  hjk>(z)  are  printed,  as  in  examples  3  and  4,  it  is  hjlHz)  that  is  used  with  j„(z)  to  determine 
the  Wronskian  check. 


1  M.  Abramowiiz  and  I  A.  Siegun,  eds.,  1965,  Handbook  of  Maihtmatkal  Functions,  U  S.  Department  of  Commerce. 
National  Bureau  of  Standards,  Washington,  DC  20402. 
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SUBROUTINE  CSP JYD ( AR , A I r N , S JR , S J I , SYR , SY I , SHR , SHI , NAK ) 

AS  OF  30  AUGUST  19S2 

DOUBLE-PRECISION  SPHERICAL  BESSEL  FUNCTIONS  FOR  COMPLEX  ARGUMENTS 
WRITTEN  BY  J.P. MASON,  ACOUSTICS  DIVISION,  NRL 
IMPLICIT  REAL*8  <A-H,0-Z) 

DIMENSION  SJR  < 1 ) ,SJI < 1) ,SYR< 1 ) ,SYI ( 1 ) ,SHR( 1 > ,SHI ( 1 ) 

LOGICAL  LT , LF 

DATA  LT/. TRUE. /,LF/. FALSE./ 

CALL  ERRSET <  72 , LF , LF , LF ,LT ,  ) 

CALL  ERRSET ( 73 , LF , LF , LF , LT , ) 

CALL  ERRSET ( 74 , LF , LF , LF , LT , ) 

C  WRITE! 5 , 104) 

104  FORMAT! '  SPHERICAL  BESSEL  FUNCTIONS  FOR  COMPLEX  ARO ' ) 

ZER0=0 . ODO 
0NE*1.0D0 
TW0=2. ODO 
THR££=3.0D0 
IZ=0 

DR=AR*AR-AI*AI 

DI=TWO*AR*AI 

CC=TWO 

EPS=  1 . 0D-1S 

WUNR*ONE 

WUNI =ZERO 

CALL  DVDD ! WUNR , WUN I , DR , D I , T 1 , T2 ) 

SRARG=DSGRTIAR*AR+AI*AI ) 

IF !SRARG.GT.0.5D0)G0  TO  29 
NP=N+1 

CALL  MLTD ! AR , A I , AR , A I , ZR , Z I ) 

ZR=ZR/TWO 
ZI=ZI/TWO 
FDNM=THREE 
HDN=ONE 
HDNM=ONE 
HDNI =ZERO 
DO  14  1=1, NP 
NN= 1-1 
EN=NN 

C  CALCULATE  Z**N/! 1X3X5. .. !2N+1 ) )  FOR  J 

IF ! NN-1 ) 2 ,6 , 3 
6  FNR=AR/ THREE 

FNI =AI /THREE 
GO  TO  5 

2  FNR=ONE 
FNI =ZERO 
GO  TO  5 

3  CALL  MLTD ( FNR , FN I , AR , A I , FNR , FN I ) 

FDNM=FDNM+TWO 

FNR=FNR/FDNM 
FNI =FNI /FDNM 
5  ANSR=ONE 

ANSI =ZERO 
PANSR=ONE 
PANSI =ZERO 
TRM=-ONE 
TIMrZERO 
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AB-QNE 

BA-THREE 

7  GNU = AB*  <  TWO#EN+BA ) 

ZRS--ZR/GNU 
ZIS--ZI /GNU 

CALL  HLTDITRM.TIM.ZRS.ZIS. TRM.TIM) 

ANSR=ANSF  -TRM 

ANS 1 = ANS I -T I M 

IF  ( ANSR . EG . ZERO ) GO  TO  15 

IF (ANSI .EQ. ZERO) GO  TO  16 

IF ( DABS ( ( PANSR-ANSR  ) /ANSR) . LE. EPS. AND. DABS <( PANSI-ANSI ) /ANSI ) 
1  . LE . EPS ) GO  TO  8 

GO  TO  17 

15  IF  ( DABS ( ( PANSI-ANSI ) /ANSI ) . LE .EPS ) GO  TO  8 
GO  TO  17 

16  IF ( DABS ( <  PANSR-ANSR ) /ANSR  > . LE .EPS > GO  TO  8 

17  PANSR-ANSR 
PANSI-ANSI 
AB-AB+ONE 
BA-BA+TMO 
GO  TO  7 

CALL  MLTD<FNR,FNI .ANSR. ANSI .SJRt I > ,SJI (I ) > 

CALCULATE  (-1X3X5. .. (2N-1 >> /Z**(N+1 >  FOR  Y 
IF (NN-1 >4.10.9 
4  GDR--ONE 

GDI -ZERO 

CALL  DUDDIGDR. GDI. AR.Al.HR. HI) 

GO  TO  11 

10  HDR-AR 
HD I* A I 

9  CALL  HLTD(HDR.HDI , AR. AI . HDR.HDI > 

HDNH=HDNH*HDN 

HDN»HDN+TMO 

CALL  DODD ( HDNM . HDN I . HDR . HD I . HR . H I > 

HR«-HR 
HI *-Hl 

1 1  ALSR»ONE 
ALSI»ZERO 
PALSR-ONE 
PALS I -ZERO 
TRN— ONE 
TIN-ZERO 
AC-ONE 
CA-ONE 

12  HNU-AC* ( CA-TWO*EN ) 

XRS— ZR/HNU 
XI9— ZI/HNU 

CALL  MLTD ( TRN .TIN. XRS .XIS.TRN.TIN) 

ALSR-ALSR-TRN 
ALSI-ALSI-TIN 
I F ( ALSR . EQ . ZERO  >  GO  TO  18 
I F ( ALS I . EQ . ZERO  >  GO  TO  19 

IF ( DABS( ( PALSR-ALSR  > /ALSR  > .LE . EPS .AND . DABS ( ( PALSI-ALSI ) /ALSI > 
1  . LE .EPS ) GO  TO  13 

GO  TO  20 

18  IF ( DABS ( ( PALSI-ALSI ) /ALSI > . LE . EPS ) GO  TO  13 
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GO  TO  20 

19  IF  <  DABS ( ( P ALSR- ALSR ) /ALSR ) .LE . EPS ) GO  TO  13 

20  PALSR=ALSR 
PALSI =ALSI 
AC=AC+ONE 
CA=CA+TWO 
GO  TO  12 

13  CALL  MLTD  <  HR  f  HI » ALSR » ALS IrSYR(I)rSYI<I)) 

C  WRITE(5,200) I ,SYR< I) ,SYI < I) 

C200  FORMAT  < 15  »2 ( 1 X  >  D23 . 16) ) 

C  IF  NAK  =2 »  PUT  Y'S  IN  SYR  AND  SYI. 

C  IF  NAK  =  3 »  STORE  Y'S!  PUT  H'S  INTO  SYR  AND  SYI. 

C  IF  NAK =5  r  PUT  Y'S  INTO  SYR  AND  SYIJ  PUT  H'S  INTO  SHR  AND  SHI. 

I F ( NAK . EQ . 2 ) GO  TO  14 
IF (NAK .EQ.5 )G0  TO  50 
YRR=SYR< I ) 

YI I =SY I < I ) 

IF(AI .LT. ZERO) GO  TO  51 
SYR (I)=SJR(I)-YII 
SYI(I)=SJI(I )+YRR 
GO  TO  14 

51  SYR  (I)«SJR(I>+YII 

SYI<I)»SJI<I)-YRR 
GO  TO  14 

50  IFIAI.LT. ZERO ) GO  TO  48 

SHR( I ) =SJR( I ) -SYI ( I ) 

SHI ( I ) =SJI ( I )+SYR( I ) 

GO  TO  14 

48  SHR( I ) =SJR< I )+SYI <  I  ) 

SHI(I)*SJI(I)-SYR<I> 

14  CONTINUE 
RETURN 

29  DSN=DSIN(AR> 

DCS«DCOS<AR) 

EXYL“D£XP ( A I ) 

EXYS*DEXP ( -A I ) 

XSN«AR*DSN 

XCO*AR*DCS 

YSN=AI*DSN 

YCO=AI*DCS 

ZXY* AR*AR+A I *A I 

TZXY*TMO*ZXY 

S JZRL*  < XSN+YCO ) /TZXY 

S JZRS* ( XSN-YCO ) /TZXY 

SJZIL* ( XCO-YSN ) /TZXY 

SJZIS* ( -XCO-YSN ) /TZXY 

SYZRL*EXYL* ( -S JZIL ) 

SYZRS=EXYS*SJZIS 
SYZIL*EXYL*S JZRL 
SYZIS*EXYS* ( -S JZRS ) 

SJZRL*EXYL*S JZRL 
SJZRS*EXYS*S JZRS 
S JZIL “EX YL*S JZIL 
SJZIS*EXYS*S JZIS 
SJR< 1 ) “SJZRL+SJZRS 
SJI  <  l  )«SJZIL-*-SJZIS 
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S JR  <  Z  >  =  <  < AR*S  JZRL+A  I-*S  JZ IL  )  /ZXY+SYZRL  )  + 

1  ( < AR*SJZRS+AI*SJZIS>/ZXY+SYZRS> 

SJI <2>  =  < < -A  I *S JZRL+AR*S JZ I L ) /ZXY+SYZIL )  + 

1  ( ( -A I#S JZRS+AR*S JZ IS ) /ZXY+SYZIS ) 

NH0  =  0 

IF(DABS(AI  )  .LT.5.0DOGQ  TO  43 

C  CALCULATE  HANKEL  FUNCTIONS  AND  THEN  THE  NEUMANN  FUNCTIONS. 

N  HO=  1 

YEX  =  DEXP ( -DABS  <  A I  )  ) 

ANUR=YEX#DSN 

ANUI = YEX*DCS 

IF ( AI . OE . ZERO  >  ANUI =-ANUI 

CALL  DODD ( ANUR , ANU I , AR , A I  . HRZ . H I Z  ) 

CALL  MLTD ( AR , AI , AR . A I , ZSR  >  ZS I > 

CALL  DODD  <  ANUR  >  ANU  I  >  ZSR  .  ZSI  »  HRW  .  HI W > 

IF ( A I ) 38 . 39 . 39 

38  ANUR=-ANUR 
GO  TO  40 

39  ANUI =-ANUI 

40  CALL  DODD ( ANU I , ANUR .AR.AI > HOA . HOB ) 

C  IF  NAK  =2 .  STORE  H'SJ  PUT  Y'S  INTO  SYR  AND  SYI. 

C  IF  NAK=3»  PUT  H'S  INTO  SYR  AND  SYI 

C  IF  NAK  =  5 .  PUT  Y'S  INTO  SYR  AND  SYI!  PUT  H'S  INTO  SHR  AND  SHI. 

IF< NAK . LT . 5 ) GO  TO  54 
SHR< 1)=HRZ 
SHI ( i >=HIZ 
SHR( 2 ) =HRW-HGA 
SHI <  2 ) =HIW-HOB 
GO  TO  55 

54  IF ( NAK . EQ . 2 ) GO  TO  58 

SYR ( 1 ) =HRZ 

SYI ( 1 )=HIZ 
S YR ( 2  >  =HRW-HOA 
SYI (2) “HIW-HOB 
GO  TO  36 

56  HRM=HRW-HOA 
HIW=HIW-HOB 

SYR (  1  )=-SJI(  n+HIZ 
SYI ( 1 ) *SJR < 1 ) -HRZ 
SYR(Z) =-SJI (2J+HIW 
SYI <  2 ) =S JR  <  2 ) -HRW 
GO  TO  57 

55  SYR( 1 )*-SJI( 1 )  +SHI ( 1 ) 

SYI < 1 > *SJR< 1 >-SHR( 1 ) 

SYR ( 2  )  =  -S JI ( 2 ) +SHI <  2 ) 

SYI(2)=SJR<2) -SHR  <  2 ) 

57  I F  <  A I . GE . ZERO ) GO  TO  36 
SYR ( 1 ) =-SYR< 1 ) 

SYI ( I ) --SYI ( 1 ) 

SYR ( 2 ) =-SYR  <  2 ) 

SYI  (2)=-SYI(2) 

GO  TO  36 

C  CALCULATE  NEUMANN  FUNCTIONS  AND  THEN  THE  HANKEL  FUNCTIONS. 

43  SYR( 1 ) -SYZRL+SYZRS 

SYI < 1 > =SYZIL+SYZIS 

SYR (2 ; =< ( AR»SYZRL+AI*SYZIL) /ZXY-SJZRL )+ 
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1  ( <AR*SYZRS+AI*SYZIS) /ZXY-SJZRS) 

SYI ( 2) * ( < -AI*SYZRL+AR#SYZIL ) /ZXY-SJZIL  >  + 

1  <  <-AI*SYZRS+AR*SYZIS>/ZXY-SJZIS> 

C  IF  NAK *2.  PUT  Y'S  INTO  SYR  AND  SYI. 

C  IF  NAK *3.  STORE  Y'S;  PUT  H'S  INTO  SYR  AND  SYI. 

C  IF  NAK =5.  PUT  Y'S  INTO  SYR  AND  SYI,'  PUT  H'S  INTO  SHR  AND  SHI. 

42  I F ( NAK . EQ . 2  >  GO  TO  36 

I F  ( NAK  .  EG .  5 )  GO  TO  52 
YRZ»SYR< 1 ) 

YIZ*SYI < 1 ) 

YRW"SYR<2) 

YIW»SYI(2) 

IF ( AI . LT . ZERO > GO  TO  53 
SYR( 1)*SJR(1)-YIZ 
SYI(1)*SJI(1) +YRZ 
SYR<2) *SJR(2)~YIM 
SYI  (2)*>SJI <2>+YRW 
GO  TO  36 

53  SYR( 1 )*3JR( 1 )+YIZ 

SYI(1)-8JI(1)'YRZ 
SYR ( 2 ) *S JR( 2 ) +YIW 
S Y I ( 2 ) «S J I ( 2 ) - YRM 
GO  TO  36 

52  IF(AI .LT .ZERO I GO  TO  41 

SHR ( 1 ) *S JR ( 1 ) -SY 1(1) 

SHI < 1 ) *8 J 1(1) +8YR  < 1 ) 

SHR ( 2 ) "S JR ( Z ) -8Y I <  2 ) 

SHI (2) *S J1 ( 2 ) ♦SYR( 2 ) 

GO  TO  36 

41  SHR ( 1 ) *8 JR  < 1 ) +SY 1(1) 

SHI ( 1 ) "SJI  ( i )-SYR( 1 ) 

SHR(2)«SJR(2)+SYI<2) 

SHI (2) *SJI (2)-SYR<2> 

36  IF ( N. LE. 1 ) RETURN 

C  THE  J's’.  Y'S.  AND  H'S  FOR  N*0  AND  N~1  HAVE  BEEN  GENERATED. 

M«N+l 

C  FIND  REMAINING  J'S. 

NN-SRARG+30 

IF ( (N+30) ,GT.NN)NN*N+30 
GDR-SJR(Z) 

GDI«SJI<2) 

30  SJR(NN)*ZERO 

SJI <NN) *ZERO 
SJR(NN-1 ) *1 .0D-20 
SJI<NN-l)«1.0D-20 
NM*NN-2 
DO  31  K  *  2 .  NM 
KK  *NN-K 

CALL  DVDD  ( S  JR  (  KK  + 1  )  .  S  J  I  <  KK  + 1  >  .  j-  1 i  KK  >  ,  S  J I  (KK  )  ) 

CALL  ERRSET <  72 .LT » LF »LF « LF . > 

CALL  ERRSNS 

SJR(KK>»<CC*KK+0NE>*SJR<KK)-SJR(KK+2> 

CALL  ERRSNS < NUN, ,. > 

I F ( NUN . EQ . 72  >  GO  TO  24 
CALL  ERRSNS 

SJI  (KK)  =  (CC*KK«-ONE>*SJI  (KK)-SJI  (KK+2) 
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CALL  ERRSET (72>LF>LF,LF.LT.  ) 

CALL  ERRSNS < NUM  .  ,  ,  ) 

IF(NUM.EQ.72)G0  TO  24 

31  CONTINUE 

CALL  DODD  C GDR  > GD I , S JR ( 2 ) . S J I < 2 )  -RAR. RAI  ) 

C  IF  THERE  WAS  AN  UNDERFLOW  IN  THE  DODD  SUBROUTINE  AND  EITHER  RAR 

C  OR  RAI  WAS  MADE  EQUAL  TO  ZERO  .  NN  SHOULD  BE  REDUCED. 

IF ( RAR . NE . ZERO . AND . RAI . NE . ZERO ) GO  TO  G7 
IF ( DABS (SJR(2> ) .LT. DABS ( S J I ( 2 ) ) ) GO  TO  68 
IF ( RAR. NE. ZERO) GO  TO  69 

IF ( GDR . EQ . ZERO .AND.SJI (2) .EQ. ZERO) GO  TO  69 
GO  TO  24 

69  I F ( RAI . NE . ZERO ) GO  TO  67 

IF(GDI .EQ. ZERO. AND.SJI (2) .EQ.ZERO)GO  TO  67 
GO  TO  24 

68  IF ( RAR. NE. ZERO) GO  TO  70 

IF ( GD I . EQ . ZERO . AND . S JR  <  2 ) . EQ . ZERO ) GO  TO  70 
GO  TO  24 

70  IF (RAI .NE. ZERO) GO  TO  67 

IF(DDR. EQ. ZERO. AND. SJR(2) ,EQ.ZERO)GO  TO  67 
GO  TO  24 

67  DO  32  K  =3  i  M 

TR=S JR ( K ) 

T I =S J I ( K  ) 

32  CALL  MLTD(TR,TI , RAR, RAI ,SJR(K>  ,SJI(K>  ) 

S JR ( 2 ) =GDR 

SJI (2) =GDI 

C  FIND  REMAINING  Y'S  AND  H'S. 

I F  <  NHO . EQ . 1 )G0  TO  44 

C  IF  NAK  =2 »  PUT  Y'S  INTO  SYR  AND  SYI. 

C  IF  NAK  =3 »  STORE  Y'SJ  PUT  H'S  INTO  SYR  AND  SYI. 

C  IF  NAK  =  5 »  PUT  Y'S  INTO  SYR  AND  SYI I  PUT  H'S  INTO  SHR  AND  SHI. 

I F ( NAK . EQ . 3 ) GO  TO  66 

22  DO  23  K=3.M 

CALL  DODD ( SYR ( K- 1  )  , SY I < K- 1 > , AR , A I ,SYR(K)  ,SYI  (K  >  ) 

SYR ( K ) = ( CC*K- THREE )*SYR(K>-SYR(K-2> 

SYI (K )= (CC*K-THREE)*SYI (K )-SYI (K-Z) 

I F ( NAK . EQ . 2 ) GO  TO  23 
I F (AI .LT. ZERO) GO  TO  45 
47  SHR(K ) =SJR(K )-SYI <K ) 

SHI (K ) =  SJI <K )+SYR(K  > 

GO  TO  23 

45  SHR (K)=SJR(K)+SYI(K) 

SHI (K ) =SJI (K)-SYR(K ) 

23  CONTINUE 
RETURN 

6G  DO  60  K=3,h 

CALL  DODD ( YRW , Y I W , AR , A I , YRT , Y I T ) 

YRT=  <CC*K -THREE >*YRT-YRZ 
Y I T  = ( CC*K- THREE )#YIT-YI2 
IF  <  AI . LT . ZERO ) GO  TO  58 
SYR ( K ) =  S JR  <  K ) -Y I T 
SYI (K  >  =SJI t  K ) + YRT 
GO  TO  53 

58  SYR ( K ) - S  JR ( K ) +Y  I  T 

SYI (K ) -SJI (K ) -YRT 
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59  YRZ=YRW 
Y I Z  =  Y I W 
YRW= YRT 
Y I W= Y I T 

60  CONTINUE 
RETURN 

C  IF  NAK=2 ,  STORE  H'S?  PUT  Y'S  INTO  SYR  AND  SYI. 

C  IF  NAK  =  3 1  PUT  H'S  INTO  SYR  AND  SYI. 

C  IF  NAK  =5  r  PUT  Y'S  INTO  SYR  AND  SYI  ;  PUT  H'S  INTO  SHR  AND  SHI. 

44  IF ( NAK . NE . 5 ) GO  TO  61 

DO  46  K*3,M 

CALL  DVDD(SHR(K-1 ) ,SHI (K-l ) , AR , AI , SHR( K ) , SHI <K) ) 

SHR ( K  )  =  ( CC#K-THREE ) «SHR ( K ) -SHR ( K-2  > 

SHI  (K )  «=  (CC#K -THREE  )  *SHI  (K)-SHi (K-2) 

SYR(K )*-SJI (K) ♦SHI ( K ) 

SYI(K)*SJR(K) -SHR  <K ) 

IF ( AI .GE . ZERO ) GO  TO  46 
SYR(K)=-SYR(K) 

SYI <  K ) =-SYI ( K ) 

46  CONTINUE 

RETURN 

61  IF ( NAK . EG . 3 ) GO  TO  62 
DO  63  K  =  3,M 

CALL  DVDD<HRW.HIW,AR, AI ,HRT,HJT) 

HRT  = ( CC*K-THREE ) *HRT-HR2 
HIT  =  <CC*K-THREE)*HIT-HIZ 
SYR(K)=-SJI(K)+HIT 
SYI(K)=SJR(K) —HRT 
I F ( A I . GE . ZERO ) GO  TO  64 
SYR ( K ) =-SYR ( K ) 

SYI (K ) =-SYI (K ) 

64  HRZ=HRW 
HIZ=HIW 
HRU=HRT 
H I W=H I T 

63  CONTINUE 

RETURN 

62  DO  65  K  =  3 . M 

CALL  DODD (SYR (K-l) . SY I ( K- 1 ) , AR, AI r  SYR ( K  > , SY I ( K ) ) 

SYR ( K ) = ( CC*K -THREE )*SYR(K>-SYR(K-Z) 

SYI ( K ) = ( CC*K- THREE )*SYI(K)-SYI (K-Z) 

65  CONTINUE 
RETURN 

24  NN=NN- 1 

WRITE ( 5 , 26 ) NN 

26  FORMAT  <1X,'  NN  REDUCED  TO  ',16) 

IZ=IZ+1 

I F  < IZ.GT. 25) RETURN 

GO  TO  30 

END 

% 
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SUBROUT  I NE  DVDD ( XA , YA , XB , YB , XC , YC  ) 

AS  OF  11  JANUARY  1983 

WRITTEN  BY  JANET  P.  MASON 
IMPLICIT  REALMS  (A-H,  -Z) 

LOGICAL  LT.LF 

DATA  LT/ .TRUE. / ,LF/ .FALSE. / 

2ER0=O . ODO 

IF  (  XB  .  NE  .  ZERO  .  OR  .  YB  .  NE  .  ZERO  )  GO  TO  3 
WRITE15, 100) 

WRITE1G, 100) 

100  FORMAT  ('  BOTH  REAL  AND  IMAGINARY  PARTS  OF  DENOMINATOR  ARE  ZERO') 
RETURN 

3  CALL  ERRSET ( 72 , LT, LF , LF , LF  »  > 

CALL  ERRSET ( 73 » LT  .  LF , LF , LF  , ) 

CALL  ERRSET ( 74 , LT , LF , LF , LF , ) 

DENOM =XB*XB+YB*YB 
IF ( DENOM . EG . ZERO ) GO  TO  1 
XX= ( XA*XB+YA*YB ) /DENOM 
I F ( XX . EG . ZERO ) GO  TO  1 
YC=(YA*XB-XA*YB) /DENOM 
IF ( YC . EG . ZERO ) GO  TO  1 
XC*XX 
RETURN 

1  CALL  ERRSET(72.LF,LF.LF,LT, ) 

CALL  ERRSET  <  73 . LF . LF , LF , LT , ) 

CALL  ERRSET < 74. LF. LF. LF. LT, ) 

IF <  DABS ( XB ) .LT . DABS ( YB ) ) GO  TO  2 
8  DC-YB/XB 
AC=XA/XB 
BC=YA/XB 

CALL  ERRSET ( 74, LT, LF, LF, LF, ) 

DENOM= 1 . 0D0+DC*DC 

IF  DC* DC  UNDERFLOWS,  DENOM  WILL  EQUAL  1 . ODO 
XC“ <  AC+BC*DC ) /DENOM 
YCX ( BC-AC*DC ) /DENOM 
CALL  ERRSET ( 74 , LF , LF , LF , LT , ) 

RETURN 

2  AD=XA/YB 
CD=XB/YB 
BD=YA/YB 

CALL  ERRSET  <  74 , LT , LF , LF , LF , > 

DENOM= 1 . ODO+CD*CD 

IF  CD*CD  UNDERFLOWS,  DENOM  WILL  EQUAL  l.ODO 
XC=(BD+AD*CD) /DENOM 
YC=(-AD+BD*CD) /DENOM 
CALL  ERRSET ( 74 , LF , LF , LF , LT , ) 

RETURN 

END 


SUBROUTINE  MLTD ( XA , YA , XB , YB , XC , YC ) 
AS  OF  31  JULY  1S7B 

WRITTEN  BY  JANET  P.  MASON 
IMPLICIT  R£AL*8  (A-H.O-Z) 
XX=XA#XB-YA#YB 
YC*XA*YB+YA*XB 
XC-XX 
RETURN 
END 
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PROGRAM  TSPHBF 

C  AS  OF  11  JANUARY  1983 

C  WRITTEN  BY  JANET  P.  MASON 

IMPLICIT  REALMS  <A-H,Q-Z> 

PARAMETER  NA=2,NB=1202 

DIMENSION  X(NA) ,Y(NA> , MAX ( NA ) , NAF ( 3 ) 

DIMENSION  A(NB) ,B(NB> ,C(NB> ,D(NB) ,E(NB) ,F(NB> 

DATA  X/IOOO. ODO, 1000. ODO/ 

DATA  Y /600 . ODO , SOO . ODO/ 

DATA  MAX/4,1167/ 

DATA  NAF/2 ,3,5/ 

ZER0=0 . ODO 

ONE=l .ODO 

EXD=1 .OD+153 

DO  5  L=3,3 

DO  1  1=1, NA 

AR=ZERO 

AI=ZERO 

NMAX=MAX ( I )  +  1 

1 F ( NAF ( L  > . NE . 5 ) GO  TO  8 

CALL  CSPJYD(X< I ) , Y ( I ) , NMAX , A , B , C , D , E , F , NAF ( L ) ) 

GO  TO  9 

8  CALL  CSPJYD(X( I) , Y< I ) , NMAX , A , B , C , D , G, H , NAF ( L ) ) 

9  CALL  MLTD ( X ( I ) , Y< I ) , X ( I ) , Y( I > , ZSR , ZSI > 

CALL  DODD (ONE, ZERO, ZSR, ZSI , ZSR, ZSI ) 

WRITE  <5,S)X(I),Y(I) , ZSR, ZSI 

G  FORMAT  (  7X  ,  'Z  =  '  ,2<  1X,D23.  IS)  , /,  '  1/(Z*Z>  =  '  ,2(  1X,D23. 16) 

1  ,//,29X, 'REAL  PART ',14X,  'IMAGINARY  PART',/) 

NNMAX=NMAX+1 
DO  4  J=1 ,NNMAX 
NN= J-2 

IF ( J .EG. 1 >GO  TO  4 
IF(NAF(L>-3)7, 10, 12 

7  IF  (  DABS  ( A  (  J  )  )  .  GT .  EXD  .  OR  .  DABS  ( B  (  J  >  )  .GT.EXD. 

1  OR.DABS(C( J-l ) ) . GT . EXD . OR . DABS ( D ( J- 1 > ) .GT.EXD. 

2  OR.DABS(A( J-l ) ) . GT .EXD .OR .DABS ( B( J-l >  > .GT.EXD. 

3  OR .  DABS  (C  (  J  )  )  .  GT  .  EXD.  OR  .DABS  ( D  (  J  >  )  .  GT  .  EXD  > GO  TO  13 
CALL  MLTD  (  A  (  J  )  ,B<  J  )  ,C<  J-l  >  ,D(  J-l  )  ,RRR,RRI> 

CALL  MLTD<A( J-l ) ,B(J-1),C(J),D(J> ,SRR,SRI > 

AR=-RRR-SRR 
AI=RRI-SRI 
C.0  TO  13 

12  BR=B  <  J ) *E  <  J- 1 >  +A  <  J ) #F  <  J- 1 ) -B  <  J- 1 >  *E ( J  > - A ( J- 1 >  *F  <  J  > 

BI=-A< J)*E< J-l )+B< J)*F( J-l )+A( J-1)*E( J)-B( J-1)*F(J) 

IF( Y( I ) .GE. ZERO) GO  TO  13 

BR=-BR 
BI=-BI 
GO  TO  13 

10  AR=B( J)*C< J-l )+A( J)*D( J-l >-B( J-l ) *C ( J ) -A (J- 1 ) *D< J ) 

AI=-A< J)#C< J-l >+B< J)*D( J-l >+A< J-l ) *C  < J)-B< J-l )*D( J) 

IF ( Y ( I ) . GE . ZERO ) GO  TO  13 

AR=-AR 
AI =-AI 

13  IF(MAX( I ) .GT.ZO.AND.NN.LT. MAX ( I >-4)G0  TO  4 
IF(NAF(L)-3) 14, 16, 18 

14  WRITE <5, 15)NN,A< J-l ) ,B< J-l ) ,C< J-l ) ,D( J-l ) ,AR,AI 
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i 


15 


16 

17 


18 

IS 


4 

2 

1 

5 


* 


1 

2 


1 

2 


1 

2 

3 


FORMAT! '  N  =  ' , I4,2X, 'SPHJ<Z>  =  '  ,2<  IX, 023. 16) ,/, 

11X,'SPHY(Z>  =  '  ,2< 1X.D23. 16) ,/, 

9X . ' WRONSK I AN  =  ' ,2 < IX ,023. 16) ) 

GO  TO  4 

WR I TE ( 5 , 17) NN  » A ( J- 1 ) , B ( J-l >  , C  <  J-l ) >D( J-l ) , AR , A I 
FORMAT ( '  N  =  '  ,14, 2X,  'SPHJ(Z)  =  ' , 2< IX , 023. 16 > , / . 

1 IX i  '  SPHH ( Z )  =  '  r2< IX, 023. 16) , / , 

9X , 'WRONSK IAN  =  ' ,2( 1X.D23. 16) ) 

GO  TO  4 

WRITE<5, 19)NN,A( J-l ) ,B< J-l ) ,C< J-l ) ,D< J-l ) ,E< J-l ) ,F( J-l ) ,BR»BI 
FORMAT (  '  N  =  ' , 14, 2X.  'SPHJ(Z)  =  ' , 2 ( 1 X , D23 . 16 ) , / , 
llX.'SPHY(Z)  =  '  ,2< IX, 023. 16) , / , 

HX.'SPHH(Z)  =  ' » 2  ( 1 X , 023 .16),/, 

9X,  'WRONSK IAN  =  '  ,2( IX, 023. 16) > 

CONTINUE 
WRITE <5, 2) 

FORMAT ( /// ) 

CONTINUE 

CONTINUE 

STOP 

END 


Example  1 


Z  « 

1/<Z*Z)  = 


-0. 1000000000000000D-02  -0. 1000000000000000D-03 
0 . 9704930889 1 285 1 6D+06  -0 . 1 9605920988 1 38420+06 


REAL  PART  IMAGINARY  PART 


N 


N 


N 


N 


0  SPHJ(Z) 
SPHY(Z) 
WRONSK IAN 

1  SPHJ(Z) 
SPHY  <  Z ) 

WRONSK I AN 

2  SPHJ(Z) 
SPHY(Z) 

WRONSK I AN 

3  SPHJ(Z) 
SPHY(Z) 

WRONSK I AN 


0 . 9999998350000079D+00 
0 . 99009850990 103060+03 
0 . 9704930889 1 285 1 90+06 
-0. 333333301000001 10-03 
-0. 97049358891 272B1D+06 
0. 970493088912B518D+06 
0. 6599999552333345D-07 
0.282441 7825519075D+10 
0.97049308891285130+06 
-0.9238094761640223D-1 1 
-0. 1 355 1 265783464 19D+ 14 
0.97049308891285160+06 


-0 . 3333333003333344D-07 
-0 . 9900995099008657D+02 
-0. 1960592098813842D+06 
-0 . 3333332336666725D-04 
0. 1960592098814092D+06 
-O. 1960592098813842D+06 
0.13333331447619120-07 
-0.87061 941 2 1960350D+09 
-0.19605920988138410+06 
-0.284761B7B8354505D-1 1 
0.57082235403167410+13 
-0. 1960592098B13843D+06 


Z 

1/(Z*Z) 


Example  2 

-0. 10000000000000000-02  -0. 1000000000000000D-03 
0 . 9704930889 12851 6D+06  -0. 1 96059209631 3842D+06 

REAL  PART  IMAGINARY  PART 


N  = 

0  SPHJ(Z) 

= 

0 

SPHH ( 2 ) 

= 

-0 

WRONSK IAN 

= 

0 

N  = 

1  SPHJ(Z) 

-0 

SPHH(Z) 

= 

0 

WRONSK IAN 

S 

0 

N  = 

2  SPHJ(Z) 

= 

0 

SPHH(Z) 

= 

-0 

WRONSK IAN 

= 

0 

N  = 

3  3PH J ( Z ) 

= 

-0 

SPHH(Z) 

2 

0 

WRONSK IAN 

2 

0 

•  SJjgawoJv/vovvv  t  •  vv 

.98009951 15508655D+02 
.97049308891 285 19D+06 
.333333301000001  ID-03 
.  19605920954807580+06 
.97049308891285180+06 
.  659999955 2333345D-07 
.  8706194 121960349D+09 
. 9704930889 1Z85 140+06 
. 923809476 1 6402 2 3D- 1 1 
.57082235403167410+13 
.  9704930889 1 285 160+06 


-0 . 3333333003333344D-07 
-0 . 9900985099343640D+03 
-0. 1 96059209881 3042D+O6 
-0. 3333332336666725D-04 
0 . 9704935BBB793949D+06 
-0. 196059209B813842D+06 
0.13333331447619120-07 
-0 .28244 178255 19075D+ 10 
-0.1 9605920988 1 384 1 D+06 
-0 . 284761 8788354505D- 1 1 
0. 1355126578346419D+14 
-0. 196059209881 3843D+06 
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Example  3 

2  *  0. 10000000000000000+04  0. 6000000000000000D+03 

1/<Z#Z)  =  0.34602076124567470-06  -0.6487889273356401D-06 


N  =  0  SPHJ(Z)  » 

SPHY(Z)  = 
SPNH(Z)  = 
WRONSK IAN  = 
N  =  1  SPHJ(Z)  * 

SPHY  <  Z )  = 

SPHH(Z)  = 
WRONSK IAN  = 
N  =  2  SPHJ(Z)  = 

SPHY(Z)  = 
SPHH(Z)  = 
WRONSK I AN  = 
N  =  3  SPHJ(Z)  = 

SPHY  <  Z  >  = 

SPHH(Z)  = 
WRONSK I AN  = 
N  =  4  SPHJ(Z)  = 

SPHY(Z)  = 
SPHH(Z)  = 
WRONSK IAN  = 


REAL  PART 

0. 1615056579356138+258 
-0.9189988821784188+256 
0 . 9538545 1 82398760-264 
0. 34602076 12456747D-06 
-0.9067180254704272+256 
-0. 1614411627841877+258 
-0 . 2063048989202653-263 
0. 34602076 12456748D-06 
-0.1613119869413142+258 
0.882186793001 1386+256 
-0.9596703805949663-264 
0 . 34602076 1 245674BD-06 
0 . 8454661476397939+256 
0.1611177608568539+258 
0 . 2064072544606 1 59-263 
0. 34602076 12456748D-06 
0.1608579340256789+258 
-0.7966475353394592+256 
0 . 9732759600894 1 94-264 
0 . 34602076 12456747D-06 


IMAGINARY  PART 

0.9189988821784188+256 
0.1615056579356138+258 
-0 . 2062840276226553-263 
-0. 64878892733564010-06 
0. 1614411627841877+258 
-0.9067180254704272+256 
-0.9557921307304426-264 
-0 . 64878892733564020-06 
-0.8821867930011396+256 
-0.1613119869413142+258 
0 . 20634624 1 72474 1 6-263 
-0 . 648788927335640 1 D-06 
-0.1611177608568539+258 
0.8454661476397939+256 
0 . 9654953095745764-264 
-0 . 648788927335640 1 D-06 
0.7966475353394592+256 
0. 1608579340256789+258 
-0 . 2064867297777066-263 
-0 . 648788927335640 1 0-06 


Example  4 

Z  =  0. 1000000000000000D+04  0 . 6000000000000000D+03 

1/(Z*Z)  *  0.34602076124567470-06  -0.64878892733564010-06 


REAL  PART  IMAGINARY  PART 


N  =  1163  SPH J  <  Z ) 
SPHY(Z) 
SPHH(Z) 
WRONSK I AN 
N  =  1164  SPHJ(Z) 
SPHY(Z) 
SPHH(Z) 
WRONSK IAN 


-0 . 3760144898599847+107 
-0.4750248660846224+107 
-0 . 4608145843562346- 1 13 
0 . 34602076 1 24567680-06 
-0.3098693455489950+107 
-0.31 33046845640833+ 1 06 
-0.8928271566952981-114 
0 . 34602076 1 2456769D-06 


0.4750248660846224+107 
-0.3760144898599847+107 
0.3825673442412044-113 
-0 . 6487889273356398D-06 
0 . 3 133046845640833+ 1 06 
-0 . 3098693455489950+107 
0.11611 06000839288- 1 1 2 
-0 . 6487889273356398D-06 


N  =  1165  SPHJ(Z) 
SPHY(Z) 
SPHH(Z) 
WRONSK I AN 
N  *  1166  SPHJ(Z) 
SPHY ( 2 ) 
SPHH(Z) 
WRONSK I AN 
N  =  1167  SPH J  <  2 ) 
SPHY  <  2 ) 
SPHH ( Z  > 
WRONSK IAN 


-0.1224447080537097+107 
0.1029806863014308+107 
0.1500954349634534-112 
0.34602076 l 2456768D-06 
-0.5900803226283226+105 
0.8191635895987654+106 
0.4407619877450161-112 
0.34602076124567680-06 
0.2800861011330237+106 
0.3146852038771679+106 
0.6270S70008025791-1 12 
0.34602078124567670-06 


-0.1029806863014308+107 
-0.1224447080537097+107 
0.1697564672546519-112 
-0 . 6487889273356397D-06 
-0.8191 635895987654+1 06 
-0 . 5900803226283226+ 1 05 
0 . 2049 1 49537982679- 1 1 3 
-0 . 64878892733563970-06 
-0.3146852038771679+106 
0.2800861011330237+106 
-0.5882652699931355-1 12 
-0 . 64878892733563980-06 
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